#include <iostream>

double epsilon=std::numeric_limits<double>::epsilon();//定义一个抽象类Function

class Function
{
public:
    //纯属虚数，计算函数值
    virtual double operator()(double x) = 0;
    virtual bool is_defined(double x){
        return true;
    }
    
    //计算导数，数值微分
    virtual double diff(double x){
        double f;
        //计算两点函数值的差
        f = (this->operator()(x) - this->operator()(x - 1e-8))/1e-8;
        return f;
             
    }
};

class EquationSolver{
public:
    //解非线性方程
    virtual double* solve()=0;
};

template <class F>
class Bisection : public EquationSolver
{
private:
    double a;
    double b;
    double delta;
    F f;
    int M;
public:
    //构造函数，接受初始区间[a, b]、最大迭代次数M、精度delta和函数f
    Bisection(double _a, double _b, int _M, double _delta, F _f)
    {
        a = _a;
        b = _b;
        delta = _delta;
        f = _f;
        M = _M;
        //函数在端点处没有根
        if (_f(_a) * _f(_b) > 0)
        {
            throw std::runtime_error("No root exists in the given interval.");
        }
        //端点处没有定义的情况
        if (!_f.is_defined(_a) || !_f.is_defined(_b))
        {
            throw std::runtime_error("Function is not defined at the endpoints.");
        }
    }
    double *solve(){
        double h = b - a;
        double u = f(a);
        double* result = new double[3];
        double c = a + h / 2;
        double w;
        int i = 0;
        for(i = 0;i < M; i++)
        {
            h = h/2;
            c = a+h;
            w = f(c);
            if(fabs(h) < delta or fabs(w) < epsilon)
            {
                break;
            }
            
            else if(w * u > 0)
            {
                a = c;
            }
        }
        result[0] = c;
        result[1] = i;
        result[2] = h;
        return result;
    }
};

template <class F>
class Newton : public EquationSolver
{
private:
    double x0;
    int M;
    F f;
public:
    // 构造函数，接受初始近似解x0、最大迭代次数M和函数f
    Newton(double _x0, int _M, F _f)
    {
        x0 = _x0;
        M = _M;
        f = _f;
       
    }

    double* solve(){
        double x = x0;
        double w;
        double *result = new double[2];
        int i = 0;
        for(i = 0;i < M; i++)
        {
            w = f(x);
            
            if(fabs(w)<epsilon)
            {
                break;
            }
            x= x - w / f.diff(x);
            
        }
        result[0] = x;
        result[1] = i;
        return result;
    }
};

// 模板类Secant，继承自EquationSolver，实现割线法求解非线性方程
template <class F>
class Secant : public EquationSolver
{
private:
  double x0;
  double x1;
  F f;
  int M;
  double delta;
public:
  // 构造函数，接受两个初始近似解x0、x1、最大迭代次数M、精度delta和函数f
  Secant(double _x0, double _x1, int _M, double _delta, F _f)
    {
      x0 = _x0;
      x1 = _x1;
      M = _M;
      delta = _delta;
      f = _f;
     
  }
  // 实现EquationSolver中的纯虚函数solve
  double* solve(){
    double x_0 = x0;
    double x_1 = x1;
    double dx;
    double w = f(x0), v;
    double s, t;
    double *result = new double[3];
    int i = 0;
    for(i = 0; i < M; i++)
    {
        v = f(x_1);
        // 判断是否满足终止条件
        if(fabs(v) > fabs(w))
        {
            t = x_1;
            x_1 = x_0;
            x_0 = t;
            v = f(x_1);
            w = f(x_0);
        }
        dx = x_1 - x_0;
        if(fabs(v) < epsilon or fabs(dx) < delta)
        {
            break;
        }
        s = dx / (v - w);
        x_0 = x_1;
        x_1 = x_1 - v * s;
        w = v;
    }
    result[0] = x_1;
    result[1] = i;
    result[2] = x_0;
    return result;
  }
};
